"""
    DataIO

A comprehensive Julia module for atmospheric data I/O operations in GSI,
ported from GSI Fortran I/O interface files for various model formats.

This module provides:
- NetCDF and HDF5 file reading and writing
- Model-specific format support (WRF, GFS, FV3, etc.)
- Background field I/O operations
- Analysis output and diagnostic file generation
- Observation data reading and processing
- Ensemble member I/O for hybrid data assimilation

# Supported Formats

## Background Fields
- **NetCDF**: Standard atmospheric model output format
- **GRIB2**: Meteorological data exchange format
- **NEMSIO**: NCEP model I/O format
- **Binary**: Native Fortran unformatted files

## Model Formats
- **WRF-ARW**: Weather Research and Forecasting model (Advanced Research WRF)
- **WRF-NMM**: WRF Nonhydrostatic Mesoscale Model
- **GFS**: Global Forecast System
- **FV3**: Finite-Volume Cubed-Sphere Dynamical Core
- **NAM**: North American Mesoscale Model

## Observation Formats
- **PrepBUFR**: Preprocessed observation format
- **BUFR**: Binary Universal Form for data Representation
- **ASCII**: Plain text observation files
- **NetCDF**: Observation diagnostic files

# Usage

```julia
using GSICoreAnalysis.DataIO

# Read background fields
bg_fields = read_background_fields("wrfout_d01_2021-01-01_00:00:00", format=:wrf_netcdf)

# Write analysis output
write_analysis_fields("analysis_2021-01-01_00.nc", analysis_state, grid)

# Read observations
obs_data = read_observations("prepbufr.2021010100", format=:prepbufr)

# Read ensemble members
ensemble = read_ensemble_members("ens_mem_????.nc", n_members=40)
```
"""
module DataIO

using LinearAlgebra
using SparseArrays
using StaticArrays
using Dates
using Printf
using NCDatasets


# Import parent types
using ..GSICoreAnalysis: AbstractAnalysisConfig
using ..GridOperations: AbstractGrid, GlobalGrid, RegionalGrid

# Export main types
export AbstractDataFormat, NetCDFFormat, GRIB2Format, NEMSIOFormat, BinaryFormat
export AbstractModelType, WRFType, GFSType, FV3Type, NAMType
export BackgroundFields, AnalysisFields, ObservationData, EnsembleData
export RadianceObservation  # Satellite radiance observation type

# Export main functions
export read_background_fields, write_analysis_fields
export read_observations, write_observations
export read_ensemble_members, write_ensemble_members
export create_diagnostic_file, write_diagnostic_data
export read_radiance_bufr, instrument_info  # Satellite radiance BUFR reading
export read_wind_bufr, WindBUFRConfig, WindObservationBundle  # Wind BUFR reading

# Data format types
"""
    AbstractDataFormat

Abstract base type for data file formats.
"""
abstract type AbstractDataFormat end

"""
    NetCDFFormat

NetCDF (Network Common Data Form) file format.
"""
struct NetCDFFormat <: AbstractDataFormat
    compression_level::Int
    chunking::Bool
    
    NetCDFFormat(; compression_level::Int=4, chunking::Bool=true) = 
        new(compression_level, chunking)
end

"""
    GRIB2Format

GRIB2 (Gridded Binary version 2) meteorological data format.
"""
struct GRIB2Format <: AbstractDataFormat
    discipline::Int
    category::Int
    parameter::Int
    
    GRIB2Format(; discipline::Int=0, category::Int=0, parameter::Int=0) = 
        new(discipline, category, parameter)
end

"""
    NEMSIOFormat  

NCEP Environmental Modeling System I/O format.
"""
struct NEMSIOFormat <: AbstractDataFormat
    endian::Symbol  # :big or :little
    
    NEMSIOFormat(; endian::Symbol=:big) = new(endian)
end

"""
    BinaryFormat

Native Fortran unformatted binary format.
"""
struct BinaryFormat <: AbstractDataFormat
    record_type::Symbol  # :sequential or :direct
    endian::Symbol       # :big or :little
    
    BinaryFormat(; record_type::Symbol=:sequential, endian::Symbol=:big) = 
        new(record_type, endian)
end

# Model types
"""
    AbstractModelType

Abstract base type for atmospheric model types.
"""
abstract type AbstractModelType end

"""
    WRFType

Weather Research and Forecasting model type.
"""
struct WRFType <: AbstractModelType
    core::Symbol  # :arw or :nmm
    version::String
    
    WRFType(; core::Symbol=:arw, version::String="4.3") = new(core, version)
end

"""
    GFSType

Global Forecast System model type.
"""
struct GFSType <: AbstractModelType
    resolution::String  # e.g., "T382", "C384"
    version::String
    
    GFSType(; resolution::String="T574", version::String="15") = new(resolution, version)
end

"""
    FV3Type

Finite-Volume Cubed-Sphere model type.
"""
struct FV3Type <: AbstractModelType
    cube_size::Int      # e.g., 384 for C384
    vertical_levels::Int
    
    FV3Type(; cube_size::Int=384, vertical_levels::Int=64) = new(cube_size, vertical_levels)
end

"""
    NAMType

North American Mesoscale model type.
"""
struct NAMType <: AbstractModelType
    grid_number::Int
    resolution_km::Float64
    
    NAMType(; grid_number::Int=218, resolution_km::Float64=12.0) = new(grid_number, resolution_km)
end

# Data structure types
"""
    BackgroundFields

Container for background (first guess) atmospheric fields.

# Fields
- `u::Array{Float64,3}`: Zonal wind velocity (m/s)
- `v::Array{Float64,3}`: Meridional wind velocity (m/s)
- `t::Array{Float64,3}`: Temperature (K)
- `q::Array{Float64,3}`: Specific humidity (kg/kg)
- `ps::Array{Float64,2}`: Surface pressure (Pa)
- `z::Array{Float64,3}`: Geopotential height (m)
- `grid::AbstractGrid`: Associated grid information
- `time::DateTime`: Valid time
- `metadata::Dict{String,Any}`: Additional metadata
"""
struct BackgroundFields
    u::Array{Float64,3}          # (nlat, nlon, nlev)
    v::Array{Float64,3}          # (nlat, nlon, nlev)
    t::Array{Float64,3}          # (nlat, nlon, nlev)  
    q::Array{Float64,3}          # (nlat, nlon, nlev)
    ps::Array{Float64,2}         # (nlat, nlon)
    z::Array{Float64,3}          # (nlat, nlon, nlev)
    grid::AbstractGrid
    time::DateTime
    metadata::Dict{String,Any}
    
    function BackgroundFields(u, v, t, q, ps, z, grid, time; metadata=Dict{String,Any}())
        # Validate dimensions
        nlat, nlon, nlev = size(u)
        @assert size(v) == (nlat, nlon, nlev) "v dimensions must match u"
        @assert size(t) == (nlat, nlon, nlev) "t dimensions must match u"
        @assert size(q) == (nlat, nlon, nlev) "q dimensions must match u"
        @assert size(ps) == (nlat, nlon) "ps dimensions must match horizontal grid"
        @assert size(z) == (nlat, nlon, nlev) "z dimensions must match u"
        
        new(u, v, t, q, ps, z, grid, time, metadata)
    end
end

"""
    AnalysisFields

Container for analysis output fields.

# Fields
- `u::Array{Float64,3}`: Analyzed zonal wind velocity (m/s)
- `v::Array{Float64,3}`: Analyzed meridional wind velocity (m/s)
- `t::Array{Float64,3}`: Analyzed temperature (K)
- `q::Array{Float64,3}`: Analyzed specific humidity (kg/kg)
- `ps::Array{Float64,2}`: Analyzed surface pressure (Pa)
- `increments::NamedTuple`: Analysis increments (analysis - background)
- `grid::AbstractGrid`: Associated grid information
- `time::DateTime`: Analysis time
- `metadata::Dict{String,Any}`: Analysis metadata (convergence info, etc.)
"""
struct AnalysisFields
    u::Array{Float64,3}
    v::Array{Float64,3}
    t::Array{Float64,3}
    q::Array{Float64,3}
    ps::Array{Float64,2}
    increments::NamedTuple{(:u, :v, :t, :q, :ps), Tuple{Array{Float64,3}, Array{Float64,3}, Array{Float64,3}, Array{Float64,3}, Array{Float64,2}}}
    grid::AbstractGrid
    time::DateTime
    metadata::Dict{String,Any}
    
    function AnalysisFields(u, v, t, q, ps, bg_fields::BackgroundFields; metadata=Dict{String,Any}())
        # Compute increments
        u_inc = u .- bg_fields.u
        v_inc = v .- bg_fields.v
        t_inc = t .- bg_fields.t
        q_inc = q .- bg_fields.q
        ps_inc = ps .- bg_fields.ps
        
        increments = (u=u_inc, v=v_inc, t=t_inc, q=q_inc, ps=ps_inc)
        
        new(u, v, t, q, ps, increments, bg_fields.grid, bg_fields.time, metadata)
    end
end

"""
    ObservationData

Container for observation data and metadata.

# Fields
- `obs_values::Vector{Float64}`: Observation values
- `obs_errors::Vector{Float64}`: Observation error standard deviations
- `obs_locations::Matrix{Float64}`: Observation locations (nobs × 3) [lon, lat, lev]
- `obs_types::Vector{Int}`: Observation type codes
- `obs_times::Vector{DateTime}`: Observation times
- `qc_flags::Vector{Int}`: Quality control flags
- `metadata::Dict{String,Any}`: Observation metadata
"""
struct ObservationData
    obs_values::Vector{Float64}
    obs_errors::Vector{Float64}
    obs_locations::Matrix{Float64}      # (nobs, 3)
    obs_types::Vector{Int}
    obs_times::Vector{DateTime}
    qc_flags::Vector{Int}
    metadata::Dict{String,Any}
    
    function ObservationData(values, errors, locations, types, times, qc_flags; metadata=Dict{String,Any}())
        nobs = length(values)
        @assert length(errors) == nobs "Error array size mismatch"
        @assert size(locations, 1) == nobs "Location array size mismatch"
        @assert size(locations, 2) == 3 "Location array must have 3 columns (lon, lat, lev)"
        @assert length(types) == nobs "Type array size mismatch"
        @assert length(times) == nobs "Time array size mismatch"
        @assert length(qc_flags) == nobs "QC flag array size mismatch"
        
        new(values, errors, locations, types, times, qc_flags, metadata)
    end
end

"""
    EnsembleData

Container for ensemble member data used in hybrid data assimilation.

# Fields  
- `members::Vector{BackgroundFields}`: Ensemble member background fields
- `mean::BackgroundFields`: Ensemble mean fields
- `perturbations::Vector{BackgroundFields}`: Ensemble perturbations (members - mean)
- `n_members::Int`: Number of ensemble members
- `metadata::Dict{String,Any}`: Ensemble metadata
"""
struct EnsembleData
    members::Vector{BackgroundFields}
    mean::BackgroundFields
    perturbations::Vector{BackgroundFields}
    n_members::Int
    metadata::Dict{String,Any}
    
    function EnsembleData(members::Vector{BackgroundFields}; metadata=Dict{String,Any}())
        n_members = length(members)
        @assert n_members > 0 "Must have at least one ensemble member"
        
        # Compute ensemble mean
        mean = compute_ensemble_mean(members)
        
        # Compute perturbations
        perturbations = [compute_perturbation(member, mean) for member in members]
        
        new(members, mean, perturbations, n_members, metadata)
    end
end

# Main I/O interface functions

"""
    read_background_fields(filename::String; format::Symbol=:netcdf, model_type::Symbol=:wrf) -> BackgroundFields

Read background atmospheric fields from file.

# Arguments
- `filename`: Path to background file
- `format`: File format (:netcdf, :grib2, :nemsio, :binary)
- `model_type`: Model type (:wrf, :gfs, :fv3, :nam)

# Returns
- BackgroundFields structure with atmospheric state variables

# Example
```julia
bg_fields = read_background_fields("wrfout_d01_2021-01-01_00:00:00", 
                                  format=:netcdf, model_type=:wrf)
```
"""
function read_background_fields(filename::String; format::Symbol=:netcdf, model_type::Symbol=:wrf)
    if format == :netcdf
        return read_netcdf_background(filename, model_type)
    elseif format == :grib2
        return read_grib2_background(filename, model_type)
    elseif format == :nemsio
        return read_nemsio_background(filename, model_type)
    elseif format == :binary
        return read_binary_background(filename, model_type)
    else
        error("Unsupported format: $format")
    end
end

"""
    write_analysis_fields(filename::String, analysis::AnalysisFields; format::Symbol=:netcdf, 
                         include_increments::Bool=true) -> Nothing

Write analysis fields to file.

# Arguments
- `filename`: Output file path
- `analysis`: AnalysisFields structure
- `format`: Output file format (:netcdf, :grib2, :nemsio, :binary)
- `include_increments`: Whether to include analysis increments in output

# Example
```julia
write_analysis_fields("analysis_2021-01-01_00.nc", analysis_fields, 
                     format=:netcdf, include_increments=true)
```
"""
function write_analysis_fields(filename::String, analysis::AnalysisFields; 
                              format::Symbol=:netcdf, include_increments::Bool=true)
    if format == :netcdf
        write_netcdf_analysis(filename, analysis, include_increments)
    elseif format == :grib2
        write_grib2_analysis(filename, analysis, include_increments)
    elseif format == :nemsio
        write_nemsio_analysis(filename, analysis, include_increments)
    elseif format == :binary
        write_binary_analysis(filename, analysis, include_increments)
    else
        error("Unsupported format: $format")
    end
end

"""
    read_observations(filename::String; format::Symbol=:prepbufr) -> ObservationData

Read observation data from file.

# Arguments
- `filename`: Path to observation file
- `format`: File format (:prepbufr, :bufr, :ascii, :netcdf)

# Returns
- ObservationData structure with observation values and metadata

# Example
```julia
obs_data = read_observations("prepbufr.2021010100", format=:prepbufr)
```
"""
function read_observations(filename::String; format::Symbol=:prepbufr)
    if format == :prepbufr
        return read_prepbufr_observations(filename)
    elseif format == :bufr
        return read_bufr_observations(filename)
    elseif format == :ascii
        return read_ascii_observations(filename)
    elseif format == :netcdf
        return read_netcdf_observations(filename)
    else
        error("Unsupported observation format: $format")
    end
end

"""
    write_observations(filename::String, obs_data::ObservationData; format::Symbol=:netcdf) -> Nothing

Write observation data to file.

# Arguments
- `filename`: Output file path
- `obs_data`: ObservationData structure
- `format`: Output format (:netcdf, :ascii)
"""
function write_observations(filename::String, obs_data::ObservationData; format::Symbol=:netcdf)
    if format == :netcdf
        write_netcdf_observations(filename, obs_data)
    elseif format == :ascii
        write_ascii_observations(filename, obs_data)
    else
        error("Unsupported observation output format: $format")
    end
end

"""
    read_ensemble_members(pattern::String; n_members::Int, format::Symbol=:netcdf, 
                         model_type::Symbol=:wrf) -> EnsembleData

Read ensemble member background fields.

# Arguments
- `pattern`: File pattern (e.g., "ens_mem_????.nc" where ???? will be replaced with member numbers)
- `n_members`: Number of ensemble members
- `format`: File format for each member
- `model_type`: Model type

# Returns
- EnsembleData structure with all ensemble members and statistics

# Example
```julia
ensemble = read_ensemble_members("ens_mem_????.nc", n_members=40)
```
"""
function read_ensemble_members(pattern::String; n_members::Int, format::Symbol=:netcdf, 
                              model_type::Symbol=:wrf)
    members = Vector{BackgroundFields}(undef, n_members)
    
    for i in 1:n_members
        # Replace pattern with member number (zero-padded to 4 digits)
        filename = replace(pattern, "????" => @sprintf("%04d", i))
        members[i] = read_background_fields(filename; format=format, model_type=model_type)
    end
    
    return EnsembleData(members)
end

"""
    write_ensemble_members(pattern::String, ensemble::EnsembleData; format::Symbol=:netcdf) -> Nothing

Write ensemble member data to files.

# Arguments
- `pattern`: Output file pattern
- `ensemble`: EnsembleData structure
- `format`: Output file format
"""
function write_ensemble_members(pattern::String, ensemble::EnsembleData; format::Symbol=:netcdf)
    for (i, member) in enumerate(ensemble.members)
        filename = replace(pattern, "????" => @sprintf("%04d", i))
        # Convert BackgroundFields to AnalysisFields for output
        # (In practice, you might want separate functions for this)
        analysis = AnalysisFields(member.u, member.v, member.t, member.q, member.ps, member)
        write_analysis_fields(filename, analysis; format=format)
    end
end

"""
    create_diagnostic_file(filename::String, grid::AbstractGrid, time::DateTime; 
                          format::Symbol=:netcdf) -> String

Create diagnostic file for GSI analysis output.

# Arguments
- `filename`: Diagnostic file path
- `grid`: Grid information
- `time`: Analysis time
- `format`: File format

# Returns
- File identifier or path for subsequent writes
"""
function create_diagnostic_file(filename::String, grid::AbstractGrid, time::DateTime; 
                               format::Symbol=:netcdf)
    if format == :netcdf
        return create_netcdf_diagnostic(filename, grid, time)
    else
        error("Unsupported diagnostic format: $format")
    end
end

"""
    write_diagnostic_data(file_id::String, varname::String, data::AbstractArray, 
                         metadata::Dict{String,Any}=Dict()) -> Nothing

Write diagnostic data to open diagnostic file.

# Arguments
- `file_id`: File identifier from create_diagnostic_file
- `varname`: Variable name
- `data`: Data array to write
- `metadata`: Variable metadata (units, description, etc.)
"""
function write_diagnostic_data(file_id::String, varname::String, data::AbstractArray, 
                              metadata::Dict{String,Any}=Dict())
    # Implementation would depend on the file format
    # For NetCDF, this would use NetCDF.jl functions
    @info "Writing diagnostic variable: $varname with size $(size(data))"
end

# Format-specific implementation functions

# NetCDF implementation functions
"""
    read_netcdf_background(filename::String, model_type::Symbol) -> BackgroundFields

Read background fields from NetCDF file.
"""
function read_netcdf_background(filename::String, model_type::Symbol)
    @info "Reading NetCDF background from: $filename"
    
    # This would use NetCDF.jl in a real implementation
    # For now, we'll create dummy data to demonstrate the interface
    
    if model_type == :wrf
        return read_wrf_netcdf_background(filename)
    elseif model_type == :gfs
        return read_gfs_netcdf_background(filename)
    elseif model_type == :fv3
        return read_fv3_netcdf_background(filename)
    else
        error("Unsupported model type for NetCDF: $model_type")
    end
end

"""
    read_wrf_netcdf_background(filename::String) -> BackgroundFields

Read WRF NetCDF background fields using the WRFReader module.

This function:
1. Opens NetCDF file using NCDatasets.jl
2. Reads coordinate variables (XLAT, XLONG, etc.)
3. Reads atmospheric variables (U, V, T, QVAPOR, PSFC, PH)
4. Handles WRF staggered grid (unstaggers winds)
5. Converts perturbation fields to total fields
6. Converts potential temperature to temperature
7. Creates grid object
8. Returns BackgroundFields

# Arguments
- `filename::String`: Path to WRF NetCDF file (wrfout_*, wrfinput_*)

# Returns
- `BackgroundFields`: Atmospheric state fields and metadata
"""
function read_wrf_netcdf_background(filename::String)
    # Use the WRFReader module for complete WRF file handling
    return WRFReader.read_wrf_netcdf(filename)
end

"""
    read_gfs_netcdf_background(filename::String) -> BackgroundFields

Read GFS NetCDF background fields.
"""
function read_gfs_netcdf_background(filename::String)
    # Similar to WRF but for GFS global format
    # Would handle spectral coefficients, global grid, etc.
    
    nlat, nlon, nlev = 181, 360, 64  # Global grid
    
    grid_config = GlobalGridConfig(nlon, nlat, nlev)
    grid = initialize_grid(grid_config)
    
    # Dummy fields
    u = randn(Float64, nlat, nlon, nlev) .* 15.0
    v = randn(Float64, nlat, nlon, nlev) .* 15.0
    t = 250.0 .+ randn(Float64, nlat, nlon, nlev) .* 20.0
    q = abs.(randn(Float64, nlat, nlon, nlev) .* 0.002)
    ps = 101325.0 .+ randn(Float64, nlat, nlon) .* 2000.0
    z = abs.(randn(Float64, nlat, nlon, nlev) .* 2000.0)
    
    time = DateTime(2021, 1, 1)
    
    metadata = Dict{String,Any}(
        "source_file" => filename,
        "model_type" => "GFS",
        "format" => "NetCDF", 
        "grid_type" => "latlon"
    )
    
    return BackgroundFields(u, v, t, q, ps, z, grid, time; metadata=metadata)
end

"""
    read_fv3_netcdf_background(filename::String) -> BackgroundFields

Read FV3 NetCDF background fields.
"""
function read_fv3_netcdf_background(filename::String)
    # FV3 uses cubed-sphere grid - more complex implementation needed
    # For now, treat as lat-lon for interface demonstration
    
    nlat, nlon, nlev = 192, 384, 64  # C192 grid mapped to lat-lon
    
    grid_config = GlobalGridConfig(nlon, nlat, nlev)
    grid = initialize_grid(grid_config)
    
    # Dummy fields
    u = randn(Float64, nlat, nlon, nlev) .* 12.0
    v = randn(Float64, nlat, nlon, nlev) .* 12.0
    t = 260.0 .+ randn(Float64, nlat, nlon, nlev) .* 15.0
    q = abs.(randn(Float64, nlat, nlon, nlev) .* 0.0015)
    ps = 101325.0 .+ randn(Float64, nlat, nlon) .* 1500.0
    z = abs.(randn(Float64, nlat, nlon, nlev) .* 1500.0)
    
    time = DateTime(2021, 1, 1)
    
    metadata = Dict{String,Any}(
        "source_file" => filename,
        "model_type" => "FV3",
        "format" => "NetCDF",
        "grid_type" => "cubed_sphere_mapped"
    )
    
    return BackgroundFields(u, v, t, q, ps, z, grid, time; metadata=metadata)
end

# Write functions (placeholder implementations)
"""
    write_netcdf_analysis(filename::String, analysis::AnalysisFields, include_increments::Bool) -> Nothing

Write analysis fields to NetCDF file.
"""
function write_netcdf_analysis(filename::String, analysis::AnalysisFields, include_increments::Bool)
    @info "Writing NetCDF analysis to: $filename"
    
    # In a real implementation, this would:
    # 1. Create NetCDF file with appropriate dimensions
    # 2. Write coordinate variables
    # 3. Write analysis fields (u, v, t, q, ps)
    # 4. Optionally write increments if include_increments=true
    # 5. Write metadata and attributes
    # 6. Close file
end

"""
    read_prepbufr_observations(filename::String) -> ObservationData

Read observations from PrepBUFR format file.

# Arguments
- `filename::String`: Path to PrepBUFR file

# Returns
- `ObservationData`: Structure containing observation values, errors, locations, types, and QC flags
"""
function read_prepbufr_observations(filename::String; config=nothing)
    @info "Reading PrepBUFR observations from: $filename"

    cfg = if config === nothing
        BUFRReader.PrepBUFRConfig()
    elseif config isa BUFRReader.PrepBUFRConfig
        config
    else
        error("config must be a BUFRReader.PrepBUFRConfig, got $(typeof(config))")
    end

    obs_data = BUFRReader.read_prepbufr(filename; config=cfg)

    if isempty(obs_data.obs_values)
        error("No observations decoded from PrepBUFR file: $filename")
    end

    return obs_data
end

"""
    write_netcdf_observations(filename::String, obs_data::ObservationData) -> Nothing

Write observation data to NetCDF file.
"""
function write_netcdf_observations(filename::String, obs_data::ObservationData)
    @info "Writing NetCDF observations to: $filename ($(length(obs_data.obs_values)) observations)"
end

# Utility functions
"""
    compute_ensemble_mean(members::Vector{BackgroundFields}) -> BackgroundFields

Compute ensemble mean from ensemble members.
"""
function compute_ensemble_mean(members::Vector{BackgroundFields})
    n_members = length(members)
    @assert n_members > 0 "Must have at least one ensemble member"
    
    # Use first member as template
    template = members[1]
    nlat, nlon, nlev = size(template.u)
    
    # Initialize mean fields
    u_mean = zeros(Float64, nlat, nlon, nlev)
    v_mean = zeros(Float64, nlat, nlon, nlev)
    t_mean = zeros(Float64, nlat, nlon, nlev)
    q_mean = zeros(Float64, nlat, nlon, nlev)
    ps_mean = zeros(Float64, nlat, nlon)
    z_mean = zeros(Float64, nlat, nlon, nlev)
    
    # Accumulate ensemble members
    for member in members
        u_mean .+= member.u
        v_mean .+= member.v
        t_mean .+= member.t
        q_mean .+= member.q
        ps_mean .+= member.ps
        z_mean .+= member.z
    end
    
    # Compute means
    u_mean ./= n_members
    v_mean ./= n_members
    t_mean ./= n_members
    q_mean ./= n_members
    ps_mean ./= n_members
    z_mean ./= n_members
    
    # Create metadata for ensemble mean
    metadata = Dict{String,Any}(
        "field_type" => "ensemble_mean",
        "n_members" => n_members,
        "source_files" => [member.metadata.get("source_file", "unknown") for member in members]
    )
    
    return BackgroundFields(u_mean, v_mean, t_mean, q_mean, ps_mean, z_mean,
                           template.grid, template.time; metadata=metadata)
end

"""
    compute_perturbation(member::BackgroundFields, mean::BackgroundFields) -> BackgroundFields

Compute ensemble perturbation (member - mean).
"""
function compute_perturbation(member::BackgroundFields, mean::BackgroundFields)
    u_pert = member.u .- mean.u
    v_pert = member.v .- mean.v
    t_pert = member.t .- mean.t
    q_pert = member.q .- mean.q
    ps_pert = member.ps .- mean.ps
    z_pert = member.z .- mean.z
    
    metadata = Dict{String,Any}(
        "field_type" => "ensemble_perturbation",
        "source_member" => member.metadata.get("source_file", "unknown")
    )
    
    return BackgroundFields(u_pert, v_pert, t_pert, q_pert, ps_pert, z_pert,
                           member.grid, member.time; metadata=metadata)
end

# Additional format-specific stubs for completeness
"""
    read_grib2_background(filename::String, model_type::Symbol) -> BackgroundFields

Read background fields from GRIB2 file.
"""
function read_grib2_background(filename::String, model_type::Symbol)
    @info "Reading GRIB2 background from: $filename"
    # Would use GRIB.jl or similar in real implementation
    error("GRIB2 reading not yet implemented")
end

"""
    read_nemsio_background(filename::String, model_type::Symbol) -> BackgroundFields

Read background fields from NEMSIO file.
"""
function read_nemsio_background(filename::String, model_type::Symbol)
    @info "Reading NEMSIO background from: $filename"
    # Would interface with NEMSIO library
    error("NEMSIO reading not yet implemented")
end

"""
    read_binary_background(filename::String, model_type::Symbol) -> BackgroundFields

Read background fields from binary file.
"""
function read_binary_background(filename::String, model_type::Symbol)
    @info "Reading binary background from: $filename"
    # Would read Fortran unformatted files
    error("Binary reading not yet implemented")
end

"""
    write_grib2_analysis(filename::String, analysis::AnalysisFields, include_increments::Bool) -> Nothing

Write analysis fields to GRIB2 file.
"""
function write_grib2_analysis(filename::String, analysis::AnalysisFields, include_increments::Bool)
    @info "Writing GRIB2 analysis to: $filename"
    error("GRIB2 writing not yet implemented")
end

"""
    write_nemsio_analysis(filename::String, analysis::AnalysisFields, include_increments::Bool) -> Nothing

Write analysis fields to NEMSIO file.
"""
function write_nemsio_analysis(filename::String, analysis::AnalysisFields, include_increments::Bool)
    @info "Writing NEMSIO analysis to: $filename"
    error("NEMSIO writing not yet implemented")
end

"""
    write_binary_analysis(filename::String, analysis::AnalysisFields, include_increments::Bool) -> Nothing

Write analysis fields to binary file.
"""
function write_binary_analysis(filename::String, analysis::AnalysisFields, include_increments::Bool)
    @info "Writing binary analysis to: $filename"
    error("Binary writing not yet implemented")
end

"""
    read_bufr_observations(filename::String) -> ObservationData

Read observations from BUFR file.
"""
function read_bufr_observations(filename::String)
    @info "Reading BUFR observations from: $filename"
    error("BUFR observation reading not yet implemented")
end

"""
    read_ascii_observations(filename::String) -> ObservationData

Read observations from ASCII file.
"""
function read_ascii_observations(filename::String)
    @info "Reading ASCII observations from: $filename"
    error("ASCII observation reading not yet implemented")
end

"""
    read_netcdf_observations(filename::String) -> ObservationData

Read observations from NetCDF file.
"""
function read_netcdf_observations(filename::String)
    @info "Reading NetCDF observations from: $filename"
    error("NetCDF observation reading not yet implemented")
end

"""
    write_ascii_observations(filename::String, obs_data::ObservationData) -> Nothing

Write observations to ASCII file.
"""
function write_ascii_observations(filename::String, obs_data::ObservationData)
    @info "Writing ASCII observations to: $filename"
    error("ASCII observation writing not yet implemented")
end

"""
    create_netcdf_diagnostic(filename::String, grid::AbstractGrid, time::DateTime) -> String

Create NetCDF diagnostic file.
"""
function create_netcdf_diagnostic(filename::String, grid::AbstractGrid, time::DateTime)
    @info "Creating NetCDF diagnostic file: $filename"
    # Would create NetCDF file with dimensions and coordinate variables
    return filename  # Return filename as file identifier
end

# Include reader submodules after core types are defined
include("WRFReader.jl")
using .WRFReader
include("BUFRReader.jl")
using .BUFRReader
include("WindBUFRReader.jl")
using .WindBUFRReader
include("RadianceInstrumentConfig.jl")
using .RadianceInstrumentConfig
include("RadianceBUFRReader.jl")
using .RadianceBUFRReader

end # module DataIO
